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Abstract 

A fully-explicit, self-consistent algebraic expression for the Reynolds stress, which is the 
exact solution to the Reynolds stress transport equation in the ‘weak equilibrium’ limit for two- 
dimensional mean flows for all linear and some quasi-linear pressure-strain models, is derived. 
Current explicit algebraic Reynolds stress models derived by employing the ‘weak equilibrium’ 
assumption treat the production-to-dissipation ( P/e ) ratio implicitly, resulting in an effective 
viscosity that can be singular away from the equilibrium limit. In the present paper, the set of 
simultaneous algebraic Reynolds stress equations are solved in the full non-linear form and the 
eddy viscosity is found to be non-singular. Preliminary tests indicate that the model performs 
adequately, even for three dimensional mean flow cases. Due to the explicit and non-singular 
nature of the effective viscosity, this model should mitigate many of the difficulties encountered 
in computing complex turbulent flows with the algebraic Reynolds stress models. 


*This research was supported by the National Aeronautics and Space Administration under NASA Contract 
No. NAS 1-19480 while the author was in residence at the Institute for Computer Applications in Science and 
Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23681-0001. 
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1 Introduction 


Since its advent, the algebraic Reynolds stress approach introduced by Rodi [1] has been viewed 
as one of the most sophisticated closure strategy at the two-equation level of turbulence model- 
ing. In this approach, the structural form of the Reynolds stress (uiu]) is taken to be self-similar 
in space and time. That is, the anisotropy of the Reynolds stress (&,_,) defined as 


UiUj 1 

= IF ~ 3 S “- 


( 1 ) 


is taken to be a constant: 


dbij - dblj + U k ^- « 0 . 


dt dt dxk ( ^ 

In the above equations, repeated indices denote summation, Ui and u; represent the mean and 
fluctuating velocity fields and K is the turbulent kinetic energy. The weak-equi librium condition 
(equation 2) is exact for homogeneous flows at equilibrium and a reasonable approximation for 
slowly-evolving flows. However, it should be borne in mind that the model obtained invoking 
this assumption will be used in inhomogeneous flows which may be far from equilibrium. 

When the Reynolds stress transport equation is subject to this ‘weak-equilibrium’ assump- 
tion, a set of simultaneous non-linear algebraic equations is obtained. Rodi [1] proposed that this 
set of non-linear equations for Reynolds stresses be solved numerically. The iterative numerical 
solution of the set of algebraic equations can be computationally expensive, nullifying the ad- 
vantage of a two-equation model over second-order closures. By presuming that the anisotropy 
stress tensor has the form dictated by representation theory, Pope [2] obtains semi-explicit so- 
lutions for the Reynolds stresses. In order to completely close the expression for the Reynolds 
stresses, Pope’s methodology requires the numerical solution of a single non-linear equation for 
the production to dissipation ratio. 

In an effort to obtain a completely explicit expression for the Reynolds stresses from the 
set on non-linear equations, Taulbee [3] and Gatski and Speziale [4], linearize the problem by 
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treating the production of kinetic energy (which is the main source of the non-linearity) as a 
known quantity to be specified. Again, using representation theory, explicit solutions to the set 
of the now-linearized algebraic equations are obtained. This results in an effective viscosity that 
is a function of the following: 

1 't = v t{~, Si], W«j), (3) 

where, e is the dissipation rate of kinetic energy, and Sy, and Wij are the normalized strain and 
rotation rate of the fluid. This type of algebraic Reynolds stress model requires that P/e be 
specified externally: the ratio is typically set at its equilibrium value [4]. While the linearization 
of the Reynolds stress equation about the equilibrium value of (P/e) ratio is reasonable if the 
flow is near equilibrium, the resulting model can be internally inconsistent when used away from 
equilibrium. Consider the example of an arbitrary two-dimensional mean flow. Let [P/e] eq be 
the equilibrium value of the ratio used to calculate the turbulent viscosity. The production to 
dissipation ratio implied by the model can then be calculated as follows: 




where S* is the dimensional strain rate. For an arbitrary flow away from equilibrium, Sij and 
W{j can be specified without restriction: there is no guarantee that the production to dissipation 
ratio calculated from the algebraic Reynolds stress model will be even close to the equilibrium 
value assumed to calculate the turbulent viscosity. 

The current algebraic Reynolds stress models are, therefore, either self-consistent but not 
fully explicit (Rodi[l], Pope[2]) or explicit but not always self-consistent. The premise of this 
brief paper is that an algebraic Reynolds stress closure model can be of practical value if and only 
if the model expression is fully explicit, self-consistent and non-singular, and hence, computable 
in situations away from equilibrium. Towards that end, an algebraic expression for Reynolds 
stress which has all the above attributes is derived. It is also demonstrated that this expression 
is indeed the exact solution to the two-dimensional mean flow Reynolds transport equation in 


2 



the weak equilibrium limit for all linear and some quasi-linear pressure strain models. This 
model can be very useful for the calculation of complex flows, especially in situations where the 
Rodi [1] algebraic equations are still being solved numerically. 


2 The non-linear algebraic Reynolds stress equations 


The exact Reynolds stress transport equation in an arbitrary non-inertial reference frame un 
dergoing a rotation with angular velocity fi t is given by 

du'iu'j 


dt 


( Crnkj^i^k ~h ^mki^ j^k) T — Pij T 4*ij T 


1 3 9 


( 5 ) 


where e t j* is the alternating tensor. The terms, respectively, are the time rate of change, ad- 
vection, Coriolis acceleration, turbulent transport, production dissipation (e^), pressure- 

strain correlation and pressure- viscous diffusion of the Reynolds stress: 


P _ — du >_— 
ox k 


_dU t 

"dx k 


( 6 ) 


_ duj duj 
£tJ V dxkdxk 


V, 




d . . duiuj , 


The production and dissipation rate of turbulent kinetic energy are, respectively, P — \Pa, 
and e = | £a . The dissipation rate tensor can be split into its isotropic and deviatoric parts as 
follows: 


£{j — T djj. 


( 7 ) 


The transport equation for the anisotropy tensor in non-dimensional time is derived from equa- 
tion (5): 


^ 1) = ~\s x3 - (b tlt S kj + S lk b kj - h mn S mn 6 tj ) 
at so 6 

— [bik{^jk + 2e m ^*ri m ) + bjk(wik T 2e m ^fi m )] + — IT-. 


(8) 


3 



In the above equation the following normalizations have been effected using the eddy turnover 
time: 


dt * = ^dt (9) 

_ 1 KdUi dUj 
v 2e [ dxjdx i ) 

_ _ 1 K f dUi dU^ 

~ 2 £ OXj dxi ] 

n* m = 

^ = f (tij-drj). 


It is easily seen that — = —2 b mn S m n- We consider the following type of quasi-linear pressure- 
strain model (that includes all linear models): 


n?. = -(c? + c\—)bij + c 2 s tJ + c 3 (b ik s jk + b Jk s ik - h mn s mn 8 %J ) + c 4 (b tk w; k + b jk w; k ), (io) 


where the C’s are numerical constants and 


W*j - ujij + ( 11 ) 

It can be shown that most of the pressure-correlation models ([5], [6], [7]) are special cases of 
equation (10) near weak equilibrium. Substitution of equation (10) into equation (8) and invok- 
ing the weak equilibrium condition leads to the following non-linear equation for the Reynolds 
stresses: 


M(C? - 2) - 2(C 1 1 + 2 )b mn S mn ] = [C 2 - + [C z - 2](b lk S jk + b jk S ik - ^b lm S lm 8 t] ()l2) 

+[C 4 - 2}(b ik W jk + bj k W ik ). 


In the above equation W{ 3 represents the total normalized vorticity given by 


Wij — U?i \j “h 


C 4 - 4 
C 4 ~ 2 




* 

m ’ 


(13) 
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For the sake of clarity, define the following: 


= 1; = C? + 2; L2 

l ^T~ u l,s T~ 1 ' 



( 14 ) 


The non-linear algebraic Reynolds stress equation now takes the simple form 

2 

b ZJ [L | — L\b mn S mn ] = ^2 Sij + + bjkSik ~ + bjkWik). (15) 

These equations describe the fixed points of the dynamical system of equations representing the 
transport of the anisotropy of the Reynolds stresses (equation 8). 


3 Fully-explicit solution 

At this stage, the present procedure departs from those in literature (e.g., Gatski and Speziale 
[4]). Rather than treat the Reynolds stress term within the square brackets on the left hand 
side (LHS) of equation (15) implicitly as has been done in the past, this term is now retained 
in its explicit form. We now appeal to the representation theory for providing the most general 
tensorial form of the anisotropy tensor in terms of the strain and rotation rate tensors. For 
details of this, now routine, procedure the reader is referred to [2], [3] and [4]. For arbitrary- 
three-dimensional mean flows the full integrity basis is composed of ten tensors. The functional 
form is too cumbersome to be of practical value [4]. It is customary to restrict consideration 
to the more tractable case of two-dimensional mean flows and use the resultant functional form 
of the Reynolds stress model expression for three-dimensional flows also. For two-dimensional 
mean flows, the general representation of the anisotropy tensor is given by (see Gatski and 
Speziale [4] for details): 

b tj = G\Sij + G 2 (SaW kj - WikSkj) + G 3 (S ik S kj - ^5 mn 5 mn ^), (16) 

where, G\ — G 3 are unknown coefficients which are functions of the constants of the pressure- 
strain model and the invariants of the strain and rotation rate tensors. In incompressible flows, 
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these invariants are 


m = SijSij, (17) 

m = WijWij. 

The objective now is to determine the unknown coefficients by using equation (15) as the con- 
straint. 


Determination of model coefficients. For two-dimensional mean flows, using equation 
(16), it is easy to show that 

bmnSmn ~ G\Tj\. (IS) 


Substitution of equations (16) and (18) into equation (15) yields after some manipulations 


[GiSij + G 2 (S ik W kj - W ik S k j) + G s {S ik S kj - ^S mn S mn )](L? - mL\G i) (19) 

= [L>2 + ^L 3 G 3 + 2T)2L4G2]Sij + 

2L 3 Gi(Si k S k j — — S mn S mn ) — 

U Gi(S ik W kj - W ik S k j). 


Comparison of the coefficients of the tensor on either side of equation (19) leads to the 
following constraint: 

G\[L® — rj\L\G\\ = L 2 + ~L 3 G 3 + 2t^2 L 4 G 2 - (20) 

O 

The coefficients of {S ik W kj - W ik S k j ) and (S ik S kj - ^S mn S mn 8 t:j ) yield: 

C = ~ L * C = — /on 

2 Ll-rnLXG,' 3 Ll-nyL\Gi (21) 

The problem is now reduced to that of determining Gi alone. From equations (20) and (21) we 

get the following cubic equation for G\\ 

a, (Li - viL'.G,) 2 = L 2 (Ll - % L\Gt) + [ |vi(i 3 ) 2 - 2, 2 (£ 4 ) 2 ]G 1 , (22) 
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which can be rewritten in the following standard form: 


{mL\) 2 G\ - (2r h L°L\)G\ + [(L°) 2 + mL\L 2 - - V i(L 3 ) 2 + 2r }2 (L 4 ) 2 ]G i - L\L 2 = 0. (23) 


It is immediately apparent that the cubic equation degenerates into a linear equation when 
?j! = 0 or Lj = 0. For these special cases G\ is given by 

l°l 2 


G] = 


and 


Gi = 


(L\Y + 2 m (L,y 

l°l 2 


, when r)i = 0, 


(L?) 2 + 7 n L\L 2 - I m(L 3 ) 2 + 2t l2 (L 4 y 


, when L\ — 0. 


(24) 


(25) 


General solution. For the general case, the calculation of G\ is not straightforward. This 
is due to fact that the cubic equation can produce multiple real roots and the choice of the 
appropriate solution may be difficult to make. The proper choice among the possible roots of 
equation (23) is necessary. 

The solution to the cubic equation (23) can be calculated following the standard procedure 
given in most mathematical handbooks. Define the following: 


P = - 
9 = 


21? ^ = L\L 2 


(mL\) 


1\2 ’ 


m LY 

[(L°) 2 + rjiL\L 2 — -r}\{L 3 ) 2 + 2rj 2 (L 4 ) 2 ]; 


(26) 


{mL\) 


a = ($-y); b= ^(2p 2 -9pq + 27r). 


The discriminant of the cubic equation can now be calculated: 

4 27 


(27) 


If the discriminant is positive, the cubic equation (23) has one real and two complex roots. The 
choice of G\ is obvious and we pick the real root: 


Gi = + (-| + VD)' 3 + - \/I))3 for D > 0. 


(28) 
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When the discriminant is negative equation (23) has three real roots given by 


= 

(?i 2) = 

Gf = 


P r. ~ Q 

-3 +2 VT CM( 3 ) - 


( 29 ) 


p . —a .6 2?r x 


P „ ~ a ,9 4tt x 

-3 +2 VT C05( 3 + T ) - 


In the above equations 0 is given by 


cos(0) = 


- 6/2 

v / -« 3 /27 ' 


(30) 


The choice of which root to pick is now less obvious. We now need a selection criterion to 
uniquely determine Gi . 


Selection criterion. The selection of an unique value of G\ is based on the following criterion. 
Consider a calculation in which the discriminant D changes sign passing through zero. It is, 
then , reasonable to require that G\ be a continuous function of D across D = 0. This requirement 
translates into the following selection criterion: 


lim G\ — lim Gi- 

D-+ 0 + £>— ► 0 - 

When the discriminant is nearly zero, we must have 


n i & 2 — « 3 

a < 0; and — = , 

4 27 ’ 


leading to 



We now need to consider two separate cases: when 6 > 0, and 6 < 0. 


(31) 


(32) 


(33) 


Case 1: 6 > 0. From equation (28) we have, 


lim Gi = = 

D-o+ 3 1 2 



(34) 
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From equation (30) we can infer the following: 


cos(#) = 


zEL 2 

^/-« 3 /27 


(35) 


implying that 0 = n. Substitution of this into equation (29) yields the following: 



G< 2) 

gS 3) 



(36) 


Therefore, the branch of the solution that will lead to G\ being continuous function of the 
discriminant is g[ 2] . 


Case 2: b < 0* For this case we have, 

lim Gi = — ^ + 

D-> 0+ 3 

We get 0 = 0 since, 

cos(0) = — 

/-« 3 /27 

From equation (29) one obtains: 

<?!" = -§ + 2y/^cos(0) = + 2 

G< 2 > = -| + 2 S cm( *L ) = _P_M 

3 V 3 v 3 ; 3 V 3 ’ 

G< 3 < = -?+2,Hcos(^) = -P + 

1 3 V 3 v 3 ’ 3 

The branch of the solution that is now picked is Gj 1 ^. 

The choice of a unique value for Gi when multiple real roots are possible is based on the 

requirement that Gj be a continuous function of the discriminant D in the neighborhood of 




2 tfl* 


-sw* 


(37) 


1*1/2 


1 


/on\ 
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D = 0. This is an important requirement, since, in the the course of computation with the 
model, if D = 0 is encountered, there should be no discontinuity in the anisotropy tensor. It 
can be shown that the discriminant D is positive when strain rate dominates over rotation rates 
(m > 7 / 2 )- For example, in the plane strain case, D is always negative and D — 0 is never 
encountered. In this case, the use of the above selection criterion is somewhat questionable. 

Second selection criterion. For the case when D is always negative, we employ a second 
selection criterion. We will require that G\ cannot not be always positive. A positive value of 
G] would correspond to a negative value for turbulent kinetic energy production since P ~ 
—G\SijSij. In a homogeneous turbulence calculation, a negative value of production (which 
represents energy transfer from chaotic small scales to organized large scales) at all times would 
constitute a violation of the second law of thermodynamics. Furthermore, if G\ is always 
positive, the model will predict only positive values for the the turbulent shear stresses in 
the case of plane strain flow. It is well known from direct numerical simulations that the 
equilibrium value of turbulent shear stress in a plane strain case is negative. Therefore, if G\ 
is always positive, one cannot have a viable equilibrium turbulence state. In Figure 1, we plot 
the behavior of the three real roots given in equation (29) over a very wide range of pi and 772 
values for the case D < 0. It is clear from the figure that and G^ are always positive and 
G is always negative. (It turns out that when D < 0, b is always positive.) Therefore, the 
only physically viable root is G^\ 

The behavior of the multiple roots given in Figure 1 is for the Speziale et al. [7] pressure- 
strain correlation model. Similar behavior is also observed for the pressure-strain correlation 
models of Launder et al. [5] and Gibson and Launder [6] (figures not shown). 
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Model for G\. The fully explicit expression for G\ can now be summarized as follows: 

L\L 2 /[(T?) 2 + 2r) 2 (L 4 )' 2 ], for = 0; 

Z/jZ/ 2 /[(T°) 2 + t)iL]L 2 — |t? 1 (L 3 ) 2 + 2 t 7 2 (T 4 ) 2 ] for L\ = 0; 

Gi = \ ~3 + ( — 2_dr v ^) 5 + ( 2 — for D > 0; 

3 + 2 V 3 C05 (|) 5 


3 T V“2_J v ~ r V~2 ~ v-L/;- 5 , ior ^ > u; ( 4 q) 

for Z) < 0 and b < 0; 
for D < 0 and b > 0 


3 1 “V 3 ^V 3 , 

-f + 2\/?<™(f + f 

The other two coefficients G 2 and G 3 can again be easily calculated from equation (21). 

In figure 2, the coefficient G \ given by equation (40) is plotted as a function of for various 
values of t] 2 . Figures 3 and 4, provide similar plots for G 2 and G 3 . It is clear from the plots that 
these coefficients and, therefore, the effective turbulent viscosity is non-singular. 


4 Model verification 


In this Section, we first provide an explicit demonstration that the derived expression is indeed 
the exact solution to the set of non-linear algebraic equations (15). Then we compare the equi- 
librium anisotropy predicted by the new model to that calculated from the Reynolds averaged 
Navier-Stokes calculations. 


4.1 Comparison with exact solution 


The basic objective of most algebraic Reynolds stress modeling procedures is to find an explicit 
solution to the set of simultaneous non-linear equations given by equation (15). In some simple 
cases, the solution to the set of equations can be directly obtained without resorting to the 
representation theory. For homogeneous shear, one can compare the model Reynolds stresses 
to those calculated directly. Any deviation of the explicit model results from the direct solution 
even if former compares better with experimental data is undesirable, for the modeling procedure 
can claim no extra source of physics. 

Consider the case when the mean velocity gradient field is given by 


du ' - <rs f 
dT,~ S 6 " 6 ’ 2 


(41) 
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The set of non-linear equations for the Reynolds stress (15) simplifies to 


& 


bn 

&22 


_2L° fa (V>f L\L 2 
SL\ n 2(L\y l 5 J S 
I I 1° L 

-(y + ^4)(T 3 - L 4 ) - (y - T 4 )(L 3 + L 4 )] + 2S(l\) 2 

b s _i±l±_ 

L? - b l2 SL\ 

b S till 
126 L\-b u SL\' 


( 42 ) 


In figure 5, the various components of the anisotropy tensor calculated directly from solving 
the above equation are compared against those computed from the explicit algebraic Reynolds 
stress derived in the previous section. As can be seen, the results are indistinguishable for all 
values of the strain rate ratio S. 


4.2 Comparison of equilibrium anisotropy 

Any algebraic model, at the very least should be able to calculate the equilibrium state of 
anisotropy of various basic homogeneous flows. For two-dimensional homogeneous flows, this 
comparison is more a test of the pressure-strain model rather than of the algebraic stress model- 
ing methodology. Clearly, the model performance will depend upon the choice of pressure- strain 
model: we select the model of Speziale et al [6]. The algebraic Reynolds stress model (ARSM) 
is compared against the Reynolds averaged Navier-Stokes (RANS) calculations of Speziale et al 
[7], and the experimental data of Tavoularis and Corrsin [8] for the case of homogeneous shear. 
The results are given in Table 1. The model agrees very well with data. In fact, the present 
ARSM prediction is identical with the ARSM calculations of Gatski and Speziale [4]. This is 
not surprising, since, the two ARSM models should indeed be identical in the equilibrium limit. 
However, they are very different away from equilibrium, the present model being self-consistent 
and non-singular. 
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We also perform comparison between the present ARSM and RANS calculations of Speziale 
et al [7] for the case of plane-strain (two-dimensional mean flow), and axisymmetric contraction 
and expansion (three-dimensional mean flows): the equilibrium values are given in Table L As 
is to be expected, the ARSM and RANS values are quite close for the case of plane-strain. In 
this case, the difference between the two is due to the fact that the RANS calculation employs 
the full non-linear pressure strain model, whereas, only the quasi-linear pressure-strain model 
is built into the ARSM. The good agreement between the two models in the case of three- 
dimensional flows is more surprising. Recall that the ARSM procedure uses an integrity basis 
that is complete only for two-dimensional flows. In the case of axisymmetric contraction, the 
equilibrium values of anisotropy predicted by the two models are extremely close. For the case of 
axisymmetric expansion, the agreement is not as good, but quite satisfactory. (In the algebraic 
model calculations, the value of the parameter SK/e is taken from that of the full Reynolds 
stress model calculations.) 

5 Summary and Conclusions 

In this brief paper, we derive the exact solution to the Reynolds stress transport equation in 
the weak-equilibrium limit for two-dimensional mean flows for all linear and some quasi-linear 
pressure-strain correlation models. This fixed point analysis of the Reynolds stress transport 
equation produces three roots. When the discriminant of the cubic fixed point equation is 
positive, two of the roots are complex and one is real. In this case, the real root clearly is the 
physically realistic fixed point of the Reynolds stress transport equation. When the discriminant 
of the fixed point equation is negative (as it happens to be for the plane strain case), all of the 
three roots are real. Two of the real roots lead to a negative value of turbulent viscosity and, 
hence, may be unphysical. The only root that leads to a positive value of turbulent viscosity also 
produces the equilibrium values of anisotropy that are consistent with RANS calculations for 
the plane strain case. We propose this expression for the Reynolds stress as a fully explicit, self- 
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consistent algebraic model for complex flows. It is shown that the model expression captures the 
equilibrium values of anisotropy quite accurately in basic three dimensional flows also. Further 
validation and extensive testing in complex flows is currently underway. 


Acknowledgements. The author would like to thank Professor Charles Speziale for the in- 
teresting discussions on some of the aspects presented in the paper. 
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TABLE OF EQUILIBRIUM VALUES 


Type of Flow 

Method 


Equilibrium Value 




b„ 

bn 

b 2I 

b 33 

Homogeneous Shear 

ARSM 

0.204 

-0.157 

-0.149 

-0.055 


RANS 

0.218 

-0.163 

-0.146 

-0.072 


Data 

0.20 

-0.15 

-0.14 

-0.06 

Plane Strain 

ARSM 

-0.224 

0.0 

0.289 

-0.064 


RANS 

-0.21 

0.0 

0.26 

-0.05 

Axisymmetric Contraction 

ARSM 

-0.228 

0.0 

0.114 

0.114 


RANS 

-0.22 

0.0 

0.11 

0.11 

Axisymmetric Expansion 

ARSM 

0.350 

0.0 

-0.175 

-0.175 


RANS 

0.30 

0.0 

-0.15 

-0.15 
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Fig. 1(c): Behavior of G^\ 
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Fig. 2: Behavior of G v 



Vi= S ij S ij 


Figure 2: G\ as a function of iji and t) 2 
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Fig. 4: Behavior of G 3 . 



7U=S ' ..S.. 

1 1 vj i] 


Figure 4: G '3 as a function of rji and r / 2 
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Fig. 5: Model vs. Exact solution 



Normalized strain— rate S’ 


Figure 5: Comparison of model calculation and exact solution for homogeneous shear flow case. 

Lines refer to model calculation and symbols represent exact solution. — 6n, - - - 6 12 , 622, 

and 633. 
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